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Abstract  We  search  for  a  percolating,  strong  subnetwork 
of  contacts  in  a  quasi-statically  deforming,  frictional  granu¬ 
lar  material.  Of  specific  interest  in  this  study  is  that  subnet¬ 
work  which  contributes  to  the  majority  of  the  total  deviator 
stress  and  is,  or  is  on  the  edge  of  being,  isostatic.  We  argue 
that  a  subnetwork  derived  from  the  minimal  spanning  trees 
of  a  graph — optimized  to  include  as  many  elastic  contacts 
as  possible  and  which  bear  normal  contact  forces  above  a 
given  threshold  delivers  such  a  network.  Moreover  adding 
the  strong  3 -force-cycles  to  the  spanning  tree  introduces  a 
level  of  redundancy  required  to  achieve  a  network  that  is 
almost  if  not  isostatic.  Results  are  shown  for  assemblies  of 
non-uniformly  sized  circular  particles  under  biaxial  compres¬ 
sion,  in  two-dimensions:  a  discrete  element  (DEM)  simula¬ 
tion  of  monotonic  loading  under  constant  confining  pressure, 
and  cyclic  loading  of  photoelastic  disks  under  constant  vol¬ 
ume. 
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1  Introduction 

Ioannis  Vardoulakis  and  his  collaborators  brought  soil 
mechanics  to  a  level  comparable  to  other  disciplines  of  con¬ 
tinuum  mechanics  and  as  a  result  enriched  both.  His  studies 
of  shear  bands  (strain  localization),  Cosserat  theory,  and  sta¬ 
bility  brought  those  subjects  into  the  mainstream  at  a  time 
when  the  numerical  analysis  community  was  struggling  with 
the  validity  of  constitutive  theories  for  frictional  media  viz. 
a  viz.  the  mathematical  well-posedness  of  associated  initial 
and  boundary  value  problems.  The  intense  interest  seen  today 
in  micropolar  theory  is  a  direct  result  of  his  work.  AT’s  last 
face-to-face  conversation  with  Ioannis  was  on  micropolar 
constitutive  models  that  explicitly  accounted  for  force  chain 
evolution  [1],  In  this  discussion,  Ioannis  raised  a  model  that 
he  developed  in  the  late  80s  in  which  he  envisaged  the  gran¬ 
ular  medium  to  be  a  ‘two-fractions  mixture’ — comprising 
‘weak  or  frail’  and  ‘strong  or  competent’  grains  [2],  This 
study  was  inspired  by  that  paper.  Using  a  complex  networks 
approach,  we  explore  other  properties  exhibited  by  these  two 
fractions  in  connection  with  the  macroscopic  stress  and  the 
structural  mechanics  concept  of  redundancy.  As  this  study 
integrates  several  key  concepts  and  developments  in  the 
physics  and  mechanics  of  granular  systems,  we  provide  first 
an  exposition  of  these  to  put  into  context  this  effort  before 
presenting  our  findings. 

Features  that  have  both  orientation  and  spatial  extent, 
exemplified  by  the  so-called  force  chains  (the  strong  grains 
in  [2]),  dominate  the  micromechanics  of  granular  media 
yet  fall  outside  the  domain  of  traditional  ‘local’  continuum 
mechanics.  In  devising  alternative  continuum  theories  based 
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on  non-local  and  micropolar  formalisms  one  must  reckon 
with  the  general  lack  of  empirical  evidence  to  account  for 
structural  evolution,  boundary  conditions  and  geometric  con¬ 
figuration  of  the  material  domain.  Data  extracted  from  dis¬ 
crete  element  method  (DEM)  simulations  and  photoelastic 
disk  experiments  [3-8],  provide  a  wealth  of  information  to 
support  non-traditional  theories.  Details  on  force  chain  evo¬ 
lution  form  a  key  outcome  that  has  profound  implications 
for  the  broad  science  of  granular  materials  and  especially  for 
constitutive  theory. 

One  of  us,  CT,  notes  that  in  the  90s  there  was  much  interest 
in  ‘strong  force  chains’  of  particles  and  their  contribution  to 
the  macroscopic  stress  that  quantifies  the  load-bearing  capac¬ 
ity  of  the  material.  However,  the  stress  tensor  is  a  function  of 
the  contact  information,  forces  and  local  coordinates,  rather 
than  the  particles  themselves.  In  this  context,  Radjai  et  al.  [9] 
showed  that,  for  2D  systems  of  rigid  disks,  the  deviatoric 
stress  was  entirely  due  to  the  strong  subnetwork  of  con¬ 
tacts  transmitting  larger  than  average  contact  forces.  This 
was  confirmed  in  [10]  from  3D  DEM  simulations  of  ‘soft’ 
spheres,  irrespective  of  the  elastic  properties  of  the  parti¬ 
cles.  Thornton  [11]  also  showed  that  the  tangential  contact 
forces  only  provided  a  small  (typically  ~  15%)  contribu¬ 
tion  to  the  deviator  stress.  The  same  was  found  in  triaxial 
compression  tests  for  various  granular  assemblies  compris¬ 
ing  irregularly-shaped  particles  [12].  Moreover,  for  general 
states  of  stress  (eri  <73),  it  was  demonstrated  in  [13] 

that  any  deviation  from  the  isotropic  stress  state,  i.e.  the 
deviatoric  stress,  was  almost  entirely  due  to  the  strong  sub¬ 
network  of  contacts  each  transmitting  larger  than  the  global 
average  force.  A  review  of  studies  focussed  on  exploring 
the  connection  between  the  strong  subnetwork  of  contacts 
and  the  macroscopic  stress  for  various  loading  conditions  is 
given  in  [14].  Overall,  these  observations  suggest  the  possi¬ 
bility  that  the  strong  contact  network,  described  by  Radjai 
and  co-workers  as  the  “solid-like  backbone”  of  the  mate¬ 
rial  [9] — may  be  isostatic  but  is  embedded  within  the  overall 
redundantly  constrained  (hyperstatic  or  statically  indetermi¬ 
nate)  particle  system.  Indeed  man-made  structures  depend 
crucially  on  redundant  supports  to  maintain  stability  [15]. 

While  the  force  chain  network  by  themselves  might  be 
close  to  if  not  isostatic,  structural  mechanics  dictates  that 
contacts  from  laterally  confining  neighbors  are  needed  to  pro¬ 
vide  these  columnar  load-bearing  force  chains  with  the  nec¬ 
essary  redundancies  to  maintain  stability.  In  past  analyses  of 
2D  and  3D  systems,  columnar  force  chains  have  been  found 
to  consistently  reside  in  self-organized  local  contact  topolo¬ 
gies  with  a  relatively  higher  level  of  connectivity  [6-8].  This 
is  evident  in  force  chain  particles  having  a  higher  average 
number  of  contacts  as  well  as  stabilizing  3-cycles  (i.e.  con¬ 
tacts  with  neighbors  which  are  themselves  in  mutual  contact) 
when  compared  to  other  particles  in  the  system.  The  process 
of  self-organization  in  these  dense  granular  systems  seems  to 


follow  ‘rules’  resembling  those  employed  in  the  construction 
of  man-made  structures.  Specifically,  the  system  in  the  sta¬ 
ble  regime  (e.g.  during  strain-hardening  in  the  biaxial  tests 
in  [6])  evolves  to  form  a  macroscopically  redundant  struc¬ 
ture  comprising,  at  the  mesoscopic  scale,  axial  load  bearing 
column-like  force  chains  which  are  laterally  supported  by 
truss-like  3-cyles.  However,  as  shown  in  [6],  the  bulk  redun¬ 
dancy  of  the  system  degrades  with  dilatation,  following  the 
loss  of  contacts  in  the  direction  of  extension,  with  the  greatest 
rate  of  decrease  in  the  average  number  of  contacts  per  par¬ 
ticle  recorded  during  the  unstable  strain-softening  regime. 
Evidence  that  redundancy  aids  material  stability  can  be  seen 
during  the  initial  force  chain  buckling  event:  although  this 
may  involve  the  collapse  of  multiple  force  chains,  this  crit¬ 
ical  event  occurs  just  prior  to  peak  shear  stress  when  the 
material  is  globally  stable  [6].  That  the  material  holds  post- 
critical  strength  or  load-carrying  capacity  amidst  collapse  of 
some  of  its  major  load  bearing  members — is  due  to  its  capa¬ 
bility  to  redistribute  its  internal  forces  and  moments  to  other 
contacts.  This  is  often  seen  in  engineering  structures  that  are 
commonly  endowed  with  sufficient  redundancies  for  stabil¬ 
ity  and  safety  as  highlighted  in  [15]:  “For  example,  when  a 
single  column  of  a  large  frame  buckles,  the  entire  frame  need 
not  collapse,  since  the  axial  force  from  this  column  can  be 
transferred  to  the  adjacent  columns.  ” 

In  experiments  and  other  DEM  simulations,  isostatic, 
hyperstatic  and  hypostatic  regimes  have  been  observed  to 
co-exist  simultaneously  in  different  spatial  regions  [7]. 
Specifically,  in  experiments  on  photoelastic  disks  subjected 
to  multiple  cycles  of  pure  quasi-static  shear  at  constant 
volume  [5,7],  shearbands  form,  and  force  chains  develop, 
strengthen  and  buckle,  all  accompanied  by  fluctuations  of 
local  packing  densities.  During  the  start  of  this  process,  the 
system  evolves  from  a  stress-free  initial  state  to  an  interme¬ 
diate  hypostatic  regime  below  jamming,  then  to  an  isostatic 
regime  near  jamming,  and  finally  arrives  at  a  hyperstatic 
state.  Under  shear  reversal  after  reaching  some  maximum 
strain,  force  networks  change  adaptively  to  the  switch  in  the 
shear  direction:  the  original  force  network  melts  away,  in 
part  through  loss  of  contacts  in  the  direction  perpendicular 
to  the  applied  compression,  and  a  new  force  network  forms 
with  force  chains  aligned  along  the  direction  of  the  applied 
compression. 

The  nature  of  jamming  has  been  the  focus  of  several 
recent  reviews  [16-18]  with  the  key  concepts  revolving 
around  states  which  can  be  described  as  hypostatic,  iso¬ 
static,  or  hyperstatic.  Here,  the  number  of  contacts  per  par¬ 
ticle  (denoted  typically  by  Z),  plays  an  important  role  in 
distinguishing  stability.  States  with  too  few  contacts,  such 
that  there  are  so-called  floppy  modes  which  allow  deforma¬ 
tion  within  the  system  without  energy  cost,  are  hypostatic. 
States  where  there  are  no  such  floppy  modes  are  hyperstatic, 
and  isostaticity  separates  hyperstatic  and  hypostatic  states. 
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Granular  materials  near  isostatic  states  can  show  critical-like 
behavior  in  terms  of  the  particle  rearrangement  to  exter¬ 
nal  perturbations  [19],  divergence  of  a  characteristic  length 
scale  to  a  point  force  response  [20,21],  an  enhancement  of 
the  number  of  normal  modes  at  zero  frequency  [22],  and 
anomalous  behaviors  near  the  special  jamming  transition 
point,  J  [17,23].  By  applying  coarse  graining,  Blumenfeld 
has  been  able  to  demonstrate  that  in  isostatic  systems,  the 
governing  system  of  equations  on  large  scales  may  be  hyper¬ 
bolic:  this  leads  to  a  natural  association  between  force  chains 
and  the  characteristics  of  the  system  of  equations  for  such 
media  [24,25].  In  the  same  paper  [24],  Blumenfeld  also  pro¬ 
posed  that  systems  that  are  close  to  but  not  identically  iso¬ 
static,  may  be  treated  as  a  mixture  of  two  phases  with  one 
phase  being  isostatic  and  the  other  consisting  of  a  connected, 
hyperstatic,  regime.  However,  it  is  not  clear  how  those  two 
phases  should  be  identified  in  different  spatial  regions. 

In  this  study,  we  explore  the  extent  to  which  a  percolating 
strong  subnetwork  of  contacts  satisfies  the  isostatic  condi¬ 
tion  in  a  DEM  simulation  of  monotonic  biaxial  loading  and 
in  an  experimental  cyclic  shear  test.  Previous  studies  have 
assumed  a  contact  force  threshold  fc  =  F / (F)  =  1  to  iden¬ 
tify  the  ‘strong’  network,  where  (F)  is  the  global  average 
force.  Arevalo  et  al.  [26]  have  investigated  the  contact  topol¬ 
ogy  of  highly  packed,  jammed  configurations  with  respect 
to  the  fc  force  threshold  (see,  also  [27]).  Here  we  also  relax 
this  constraint  and,  in  addition,  search  for  an  isostatic  per¬ 
colating  subnetwork  that  maximizes  the  contribution  to  the 
deviator  stress.  The  algorithm  we  propose  allows  us  to  sep¬ 
arate  the  contact  networks  into  two  parts — a  strong  network 
and  a  weak  supporting  network.  To  do  so,  we  use  two  meth¬ 
ods:  one  based  on  the  value  of  R  defined  in  (2)  of  Methods 
(and  elsewhere  [7])  and  the  other  using  the  average  force  as 
the  threshold.  The  former  produces  a  subnetwork  which  is 
exactly  at  the  isostatic  state  at  R  —  1 .  However,  we  empha¬ 
size  here  that  our  method  is  general;  it  can  be  applied  to  any 
hyperstatic  system,  whether  the  system  itself  is  near  isostatic 
state  or  far  from  the  isostatic  state. 


2  Methods 

An  assembly  of  granular  particles  can  be  represented  by  a 
mathematical  graph  or  complex  network,  where  the  network 
nodes  correspond  to  particles  and  the  links  to  contacts.  The 
rheological  response  of  the  material  to  loading  is  encoded  in 
evolving  properties  of  this  contact  network  [28].  The  strong 
(weak)  contacts  of  a  network  are  those  contacts  that  carry 
a  normal  contact  force  magnitude  above  (below)  a  given 
threshold,  where  traditionally,  the  dividing  point  between 
strong  and  weak  is  taken  to  be  the  global  mean,  or  in  nor¬ 
malized  form,  fc  =  1 .  The  stress  is  given  (here  for  circular 
particles)  by 


Oij  =  - - —  X  "?CfjAc(RAc  +  RBC)  (1) 

^PsNr  V  ceNf 

where  Np  is  the  set  of  particles  in  the  network,  V  p  is  the  local 
void  volume  of  a  particle,  Nj:  represents  the  number  of  con¬ 
tacts  in  the  subnetwork  L ,  f,  is  the  contact  force  between 
particle  A  and  particle  B,nfc  is  the  unit  branch  vector  of  the 
contact  and  RAc  and  RBc  are  the  radii  of  the  contacting  par¬ 
ticles.  Equation  (1)  gives  the  stress  for  the  solid  phase  of  the 
assembly  rather  than  the  whole  volume,  the  two  differing  by 
a  multiplicative  factor  equal  to  the  solids  fraction.  The  per¬ 
colating  subnetwork  shares  the  same  physical  space  as  the 
whole  assembly,  thus  for  purposes  of  comparison,  we  sim¬ 
plify  our  computation  by  using  solid  fraction  stresses.  The 
deviator  stress  (i.e.  second  invariant  of  the  deviatoric  stress 
tensor)  is  thus  given  by  the  difference  in  the  eigenvalues  of 
the  stress  tensor,  i.e.  D  —  4[max(k,)  —  min(k,)],  where  A.,- 
are  the  eigenvalues  of  a ij . 

An  assembly  of  grains  in  equilibrium  is  isostatic  if  there 
is  exactly  enough  contacts  for  force  and  torque,  or  moment, 
balance:  i.e.  the  number  of  unknown  independent  forces  and 
moments  is  equal  to  the  number  of  equilibrium  equations. 
If  the  assembly  is  under-constrained  (over-constrained),  i.e. 
the  number  of  unknown  independent  forces  and  moments  is 
less  (greater)  than  the  number  of  equilibrium  equations,  then 
the  network  is  hypostatic  ( hyperstatic ).  The  redundancy  of  a 
granular  assembly  can  be  quantified  using  a  procedure  dis¬ 
cussed  in  detail  in  [7]  but  key  elements  are  repeated  here  for 
the  benefit  of  the  reader.  A  scalar  ratio  R  can  be  defined  where 
R  <  1  corresponds  to  a  hypostatic  system,  R  =  1  to  an  iso¬ 
static  system,  and  R  >  1  indicates  a  hyperstatic  system.  The 
definition  of  R  accounts  for  the  number  and  types/modes  of 


Dual  support  provided  by  3-cycles 

Fig.  1  An  example  of  the  dual  supporting  role  to  force  chains  that  3- 
cycles  provide  (i.e.  resistance  to  relative  rotations  at  contact  and  lateral 
support  to  ‘prop-up'  force  chains).  If  the  normal  contact  force  carried 
at  each  edge  of  the  triangle  is  above  average  (or  a  prescribed  force 
threshold)  then  the  cycle  is  a  3-force-cycle.  If  both  3-cycles  above  are 
3-force-cycles  then  the  union  of  contacts  in  each  cycle  constitutes  the 
3-force-cycle  network 
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contact  in  a  network.  The  contacts  are  classified  according  to 
whether  they  are  elastic  or  plastic  with  a  further  distinction 
accounting  for  rolling  resistance  (contact  moment)  in  assem¬ 
blies  of  circular  disks.  In  particular,  a  stick  contact  is  such 
that  both  the  tangential  force  and  contact  moment  are  elastic 
and  independent  of  the  normal  force;  by  contrast,  a  sliding 
contact  is  such  that  the  tangential  force  is  at  the  Coulomb 
threshold  and  thus  coupled  to  the  normal  force  and  a  rolling 
contact  is  such  that  the  contact  moment  is  at  its  analogous 
Coulomb  plastic  threshold  and  thus  also  coupled  to  the  nor¬ 
mal  force.  A  sliding  and  rolling  contact  has  both  tangential 
force  and  contact  moment  at  their  respective  Coulomb  plas¬ 
tic  thresholds.  Specifically,  in  2D  R ,  the  ratio  of  the  number 
of  independent  forces  and  moments  to  the  number  of  equi¬ 
librium  equations,  is  given  by 

D  3iVstick  +  2  A slide  +  2Aroll  +  Asiide+roll 

K  =  - ,  (Z) 

3  Aphides 


DEM:  Span+3FC,  fc=1 


-2  0  2 

RS:  Span+3FC,  fc=1 


where  AstiCk,  Asiide,  Aron  and  Asiide+roii  are  the  number  of 
stick,  sliding,  rolling,  and  sliding+rolling  contacts  in  the  sub¬ 
network.  Aparticies the  number  of  particles  in  the  subnet¬ 
work.  The  pre-factors  in  (2)  represent  the  number  of  degrees 
of  freedom  needed  to  define  each  type  of  contact.  For  exam¬ 
ple,  a  stick  contact  is  below  the  limiting  state  for  all  modes 
of  deformation  which  requires  three  degrees  of  freedom  in 
2D.  In  contrast  if  the  contact  is  sliding  (rolling),  a  degree 
of  freedom  is  removed  because  the  tangential  component  of 
the  contact  force  (contact  moment)  is  coupled  to  the  normal 
force,  thus  the  multiple  is  two  [7]. 

A  subnetwork  percolates  the  material  domain  if  its  constit¬ 
uent  set  of  contacts  extend  from  one  boundary  of  the  domain 
to  the  other.  A  spanning  tree  of  a  contact  network  is  an  acy¬ 
clic  subset  of  contacts  that  connects  all  the  nodes.  A  minimal 
spanning  tree  is  a  spanning  tree  with  minimum  total  path 
length  or  contacts.  If  we  assign  a  weight  to  each  contact  or 


DEM:  Span+3FC,  R=1 


-2  0  2 


RS:  Span+3FC,  R=1 


Fig.  2  Color  online.  Snapshots  of  the  percolating  minimal  spanning  trees  and  3-force-cycle  subnetworks  (Span+3FC)  for  DEM  at  axial  strain  0.04 
(upper)  and  for  cyclic  shear  system  at  strain  step  417  (lower).  Contacts  are  blue ,  3-force-cycles  contacts  are  red  triangles  and  force  chain  nodes  are 
green 
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link,  it  is  a  spanning  tree  with  minimum  sum  of  contact 
weights.  We  weight  the  contact  matrix  according  to  contact 
mode:  1  for  a  stick  contact,  2  for  a  rolling  contact,  3  for  a  slid¬ 
ing  contact,  and  4  for  the  sliding+rolling  contacts.  Thus  we 
bias  the  search  for  a  minimal  spanning  tree  towards  finding 
stick  and  rolling  contacts — the  typical  contact  modes  within 
the  force  chain  network  [1,29]. 

The  3 -force-cycles  within  a  network  are  the  ‘strong’ 
3 -cycle  motifs:  the  force  (or,  here,  its  normal  component) 
at  each  contact  of  a  3-force-cycle  has  a  magnitude  that  is  at 
least  the  prescribed  threshold.  In  [6]  the  3-force  cycles  were 
introduced  and  discovered  to  play  an  important  supporting 
role  to  force  chains  during  and  in  the  location  of  shear  band¬ 
ing.  As  illustrated  in  Fig.  1,  the  3-force  cycles  provide  dual 
resistance  to  force  chain  buckling  by:  (1)  frustrating  particle 
rotations  crucial  for  buckling  and  (2)  propping-up  the  force 
chain  particles.  The  3-force-cycle  network  is  the  collection 


0  0.05  0.1 


Stress 


of  all  such  3 -cycles  (see,  the  red  triangle  contacts  within  the 
subnetworks  of  Fig.  2  for  examples  of  such  networks). 

Our  algorithm  for  finding  a  percolating  and  isostatic  sub¬ 
network  proceeds  in  the  following  sequence: 

1.  Set  the  normal  contact  force  magnitude  threshold  (e.g. 

fc  =  1). 

2.  Prune  all  contacts  bearing  a  normal  force  magnitude  less 

than  fc. 

3.  From  the  remaining  contact  subnetworks  of  step  2 

(a)  Find  the  3 -force-cycle  network. 

(b)  Find  the  minimal  spanning  trees  with  the  contacts 
weighted  by  contact  type. 

(c)  Construct  the  subnetwork  comprising  the  union  of 
contacts  in  the  3 -force-cycle  network  and  all  the 
minimal  spanning  trees 


R=1 


Ratio 


Fig.  3  Color  online.  DEM  system — Top  left'.  The  redundancy  of  the  networks  at  fc  =  1.  Top  right:  the  force  threshold  required  such  that  the  sub¬ 
network  (Span+3FC)  is  isostatic.  Bottom  left:  the  deviator  stress  fc  =  1  and  R  =  1.  Bottom  right:  contribution  to  deviator  stress  of  the  subnetworks 
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4.  If  we  require  a  force  threshold  such  that  the  subnetwork 
is  isostatic,  i.e.  R  —  1,  update  the  force  threshold  and 
return  to  step  2.  The  threshold  can  be  updated  using  the 
Nelder-Mead  simplex  method  to  minimize  the  cost  func¬ 
tion  (R  -  l)2. 

For  the  resulting  subnetwork  of  contacts,  we  compute  the 
redundancy  and  its  contribution  to  the  total  deviator  stress. 
Almost  all  of  the  force  chain  particles,  as  determined  from 
the  algorithm  used  in  [30],  will  be  nodes  for  a  force  threshold 
of  fc  =  1  or  less.  Note  that  some  particles  may  be  missed 
as  the  force  chain  algorithm  of  [30]  considers  only  particles 
with  above-average  particle  load  vector  magnitude  and  not 
above-average  normal  contact  force.  That  is,  the  algorithm 
to  find  force  chains  within  a  granular  assembly  in  equilib¬ 
rium  identifies  groups  of  three  or  more  contacting  particles 


fc=1 


Stress 


whose  particle  load  vectors  are:  (1)  in  quasi-linear  arrange¬ 
ment  (i.e.  consecutive  vectors  align  to  within  a  prescribed 
tolerance  angle),  and  (2)  each  with  a  magnitude  that  is  above 
the  global  average  (see  [30]  for  full  details). 

It  is  very  difficult  to  find  a  force  threshold  such  that  the 
collection  of  minimal  spanning  trees  is  isostatic.  The  redun¬ 
dancy  R  from  (2)  for  fc  =  1  and  almost  all  other  values 
above  fc  =  1  are  below  1.0.  We  must  therefore  add  some 
contacts  to  the  spanning  trees  to  increase  the  redundancy. 
Recent  studies  of  the  topology  of  laterally  supporting  con¬ 
tacts  around  force  chains  from  DEM  simulations  and  physi¬ 
cal  experiment  show,  on  average,  that  force  chains  not  only 
have  a  higher  number  of  contacts  but  also  a  higher  number 
of  3-cycles  per  particle  [6,7].  Accordingly,  we  consider  the 
union  of  the  spanning  trees  and  the  subnetwork  of  3-force- 
cycles,  the  strong  3-cycles,  for  the  given  threshold. 
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Ratio 


Fig.  4  Color  online.  Cyclic  shear  system — Top  left'.  The  redundancy  of  the  networks  at  fc  =  1 .  Top  right',  the  force  threshold  required  such  that 
the  subnetwork  (Span+3FC)  is  isostatic.  Bottom  left',  the  deviator  stress  fc  =  1  and  R  =  1.  Bottom  right :  contribution  to  deviator  stress  of  the 
subnetworks 
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3  Results 

We  present  results  of  applying  the  above  algorithm  to  two 
data  sets,  both  2D  systems  of  polydisperse  frictional  circular 
particles,  discussed  extensively  elsewhere.  The  first  is  from 
a  DEM  simulation  of  biaxial  compression  with  constant  con¬ 
fining  pressure  [1,6,7,28].  The  second  is  from  experiments 
using  photoelastic  disks  subject  to  cyclic  shear  under  con¬ 
stant  volume  [5,7],  Figure  2  (top  row)  shows  the  percolat¬ 
ing  subnetworks,  i.e.  union  of  minimal  spanning  trees  and 
3-force-cycles,  for  fc  =  1  and  R  =  I  in  the  DEM;  the 
bottom  row  is  for  one  strain  step  of  the  cyclic  shear  experi¬ 
ment.  Quantities  of  interest  for  the  subnetworks  in  the  DEM 
simulation  are  shown  in  Fig.  3.  We  see  for  f<  —  I ,  the  redun¬ 
dancy  is  typically  less  than  1 .0  throughout  loading  and  that 
the  force  threshold  required  for  R  =  1  leads  to  normalized 
values  less  than  1.0.  The  proportion  of  force  chain  parti¬ 
cles  in  these  networks  at  fc  =  1  and  the  isostatic  threshold 
are  above  99%  throughout  loading  (not  shown).  Also,  the 
deviator  stresses  of  the  subnetworks  are  highly  correlated 
with,  and  are  responsible  for,  a  majority  (typically  >80%) 
of  the  total  deviator  stress.  Figure  4  shows  results  for  the 
experimental  cyclic  shear  system.  Due  to  the  larger  errors 
of  determining  contacts  and  contact  forces,  it  is  not  always 
possible  to  find  a  rich  enough  starting  network  for  all  strain 
steps,  and  some  are  left  blank.  Furthermore,  distinguishing 
stick,  sliding,  rolling  and  sliding  and  rolling  modes  of  con¬ 
tact  is  slightly  more  involved  than  checking  contact  force  and 
contact  moment  magnitudes.  For  surviving  contacts  between 
two  disks  across  a  strain  step  we  track  the  position  of  the 
contact  on  each  disk  and  the  rotation  of  each  disk.  If  either 
disk  shows  a  non-zero  rotation  the  contact  is  either  a  roll¬ 
ing  or  a  sliding  and  rolling  mode.  If  the  change  in  position 
of  the  contact  on  either  disk  is  only  due  to  relative  motion 
then  the  contact  is  rolling  otherwise  it  is  classed  as  slid¬ 
ing  and  rolling.  If  neither  disk  rotates  but  there  is  a  change 
in  position  of  contact  then  the  mode  of  contact  is  sliding. 
A  contact  is  classified  as  stick  if  none  of  the  conditions  for 
slide,  roll,  slide  and  roll  is  met  [7].  The  resulting  subnet¬ 
works  capture  about  70%  of  the  force  chain  particles  for 
fc=  1.  This  rises  to  an  average  of  92%  for  the  isostatic 
subnetworks.  The  force  threshold  needed  for  the  isostatic 
condition  of  R  =  1,  if  it  exists,  can  be  very  low,  but  can  also 
be  above  the  usual  force  threshold  fc=  1 .  Again  throughout 
the  loading,  the  deviator  stress  of  the  subnetworks  are  highly 
correlated  with  and  can  capture  60%  (for  fc  =  1)  and  80% 
(for  R  =  1)  of  the  macroscopic  deviator  stress. 

4  Conclusion 

We  have  presented  a  method  for  finding  a  percolating,  strong 
subnetwork  of  contacts  in  a  quasi-statically  deforming  granu¬ 


lar  material  that  contributes  to  the  majority  of  the  total  devia¬ 
tor  stress  and  is  isostatic  in  the  sense  of  R  defined  in  (2)  being 
equal  to  one.  By  a  process  of  optimization  that  considers  all 
contacts  bearing  a  normal  contact  force  magnitude  above  a 
given  global  force  threshold,  we  find  that  the  union  of  the 
3-force-cycles  and  all  the  minimal  spanning  trees,  optimized 
to  include  as  many  elastic  contacts  as  possible,  delivers  such 
a  subnetwork.  If  the  threshold  is  set  to  the  global  average 
normal  contact  force,  the  resulting  subnetwork  will  typically 
percolate  across  the  system,  include  the  majority  of  the  force 
chain  particles,  and  contribute  to  the  majority  of,  as  well  as 
being  highly  correlated  with,  the  total  deviator  stress.  How¬ 
ever,  this  subnetwork  is  just  below  the  isostatic  condition  of 
R  —  1 .  If  the  isostatic  condition  of  R  —  1  is  used  as  a  con¬ 
straint,  then  the  force  threshold  required  to  achieve  this  is  less 
than  fc  =  1  for  the  DEM  but  could  be  above  fc  =  1  for  stages 
of  the  deformation  in  the  experimental  system.  The  finding 
that  granular  materials  are  inherently  bimodal,  by  which  a 
percolating  isostatic  subnetwork  of  the  material  can  be  dis¬ 
tinguished  from  the  milieu  of  more  lightly  loaded  particles, 
opens  new  avenues  for  future  study.  In  particular,  this  has 
demonstrated  how  we  can  extract  new  insights  from  DEM 
and  experiments  to  directly  facilitate  the  continued  develop¬ 
ment  of  structural  mechanics  models  of  force  chain  evolution 
such  as  that  in  [1],  and  of  predictive  continuum  models  such 
as  those  proposed  for  isostatic  systems  in  [24,25]  and  for  the 
more  general  case  in  [31]. 
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